{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolative Seperable Density Fitting (ISDF) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ISDF implementation currently provides a THC-like factorization of the two electron integrals which should converge to the FFTDF representation of the ERIs in the limit of large THC rank. This differs from the assumption of using RSGDF throughout the rest of the resource estimation scripts. However, we typically are only interested in ISDF as an initial guess for the THC factors which are then subsequently reoptimized to regularize $\\lambda$. The assumption here is that FFTDF / ISDF is a good enough approximation to the RSGDF ERIs and thus serves as a good initial guess.\n",
    "\n",
    "Let's start by comparing the ISDF-MP2 energy as a function of the THC rank parameter. Recall that $M = c_\\mathrm{THC} N/2$, where $c_\\mathrm{THC}$ is the THC rank parameter and $N$ is the number of spin orbitals. $M$ is what we call num_thc in the code. \n",
    "\n",
    "It's important to recall what we are doing in the ISDF algorithm, that is we solve\n",
    "\n",
    "\\begin{equation}\n",
    "\n",
    "u_{p\\mathrm{k}}^*(\\mathbf{r}_i) u_{q\\mathbf{k}'}(\\mathbf{r}_i) = \\sum_\\mu^M \\xi_\\mu(\\mathbf{r}_i) u_{p\\mathrm{k}}^*(\\mathbf{r}_\\mu) u_{q\\mathbf{k}'}(\\mathbf{r}_\\mu)\n",
    "\n",
    "\\end{equation}\n",
    "\n",
    "for $\\xi_\\mu(\\mathbf{r}_i)$ given a set of interpolating points $(\\{r_\\mu\\})$ which are selected from the original real space $(\\{\\mathbf{r}_i\\})$ (FFT) grid of size $N_g$ using the KMeans-CVT algorithm.\n",
    "\n",
    "For the purposes of this notebook it is helpful to use a value of $N_g$ which is smaller than that required to fully converge the FFTDF error. We will investigate this more at the end of the tutorial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "from ase.build import bulk\n",
    "\n",
    "from pyscf.pbc import gto, scf\n",
    "from pyscf.pbc.tools import pyscf_ase\n",
    "from pyscf.pbc.mp import KMP2 \n",
    "\n",
    "\n",
    "ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
    "cell = gto.Cell()\n",
    "cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
    "cell.a = ase_atom.cell[:].copy()\n",
    "cell.basis = \"gth-szv\"\n",
    "cell.pseudo = \"gth-hf-rev\"\n",
    "# Use a smaller value of Ng that would otherwise be suggested (~ 26^3) for\n",
    "# expediency + to allow exact ISDF factorization for comparison.\n",
    "cell.mesh = [11]*3\n",
    "cell.verbose = 0\n",
    "cell.build()\n",
    "\n",
    "kmesh = [1, 1, 3]\n",
    "kpts = cell.make_kpts(kmesh)\n",
    "num_kpts = len(kpts)\n",
    "mf = scf.KRHF(cell, kpts)\n",
    "mf.kernel()\n",
    "print(\"SCF energy: \", mf.e_tot)\n",
    "\n",
    "# converged SCF energy with appropriate Ng = -10.388904514046914, mesh = 28^3\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's find the ISDF THC factors using the KMeans-CVT algorithm to find the interpolating points. It's easiest to use the helper function `solve_kmeans_kpisdf` which will perform the necessary steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import openfermion.resource_estimates.pbc.thc as kthc\n",
    "# Let's use the whole real space grid for some correctness checks first.\n",
    "num_thc = np.prod(cell.mesh)\n",
    "kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
    "print(kpt_thc.__dict__.keys())"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the kpt_thc class has 4 attributes, `chi`, `xi`, `zeta` and `G_mapping`. `chi` corresponds to the cell periodic part of the Bloch orbital (i.e. $u_{p\\mathbf{k}}(\\mathbf{r}_\\mu))$. `xi` corresponds to $\\xi_{\\mu}(\\mathbf{r})$ in Eq. (1) above. To understand `zeta` and `G_mapping` it is helpful to recall we want to build\n",
    "\n",
    "$$\n",
    "\n",
    "(p\\mathbf{k}pq\\mathbf{k}-\\mathbf{Q}|r\\mathbf{k}'-\\mathbf{Q} s\\mathbf{k}') = \\sum_{\\mu\\nu} u^*_{p\\mathbf{k}}(\\mathbf{r}_\\mu))u_{p\\mathbf{k}-\\mathbf{{Q}}}(\\mathbf{r}_\\mu) \\zeta_{\\mu\\nu}^{\\mathbf{Q}\\Delta \\mathbf{G}_{\\mathbf{Q}\\mathbf{k}-\\mathbf{Q}}\\Delta\\mathbf{G}_{\\mathbf{Q}\\mathbf{k}'-\\mathbf{Q}}} u^*_{p\\mathbf{k}'}(\\mathbf{r}_\\nu)u_{p\\mathbf{k}-\\mathbf{Q}}(\\mathbf{r}_\\nu)\n",
    "\n",
    "$$\n",
    "\n",
    "So `zeta` corresponds to $\\zeta_{\\mu\\nu}^{\\mathbf{Q}\\Delta \\mathbf{G}_{\\mathbf{Q}\\mathbf{k}-\\mathbf{Q}}\\Delta\\mathbf{G}_{\\mathbf{Q}\\mathbf{k}'-\\mathbf{Q}}}$ above, and `G_mapping` is a 2D array yielding the appropriate $\\Delta G$ index given an index for $\\mathbf{Q}$ and $\\mathbf{k}$.\n",
    "\n",
    "Let's look at an example to see that everything is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.resource_estimates.pbc.hamiltonian import build_momentum_transfer_mapping \n",
    "\n",
    "momentum_map = build_momentum_transfer_mapping(cell, kpts)\n",
    "num_spatial_orbs = mf.mo_coeff[0].shape[-1]\n",
    "\n",
    "# Pick a particular momentum transfer\n",
    "Q_indx = 1\n",
    "k_indx = 1\n",
    "k_prime_indx = 0\n",
    "k_minus_Q_indx = momentum_map[Q_indx, k_indx]\n",
    "k_prime_minus_Q_indx = momentum_map[Q_indx, k_prime_indx]\n",
    "\n",
    "eri_kindices = [k_indx, k_minus_Q_indx, k_prime_minus_Q_indx, k_prime_indx]\n",
    "eri_mos = [mf.mo_coeff[kindx] for kindx in eri_kindices]\n",
    "eri_kpts = [mf.kpts[kindx] for kindx in eri_kindices]\n",
    "\n",
    "eri_exact = mf.with_df.ao2mo(eri_mos, eri_kpts, compact=False).reshape((num_spatial_orbs,)*4)\n",
    "\n",
    "\n",
    "kthc_eri_helper = kthc.KPTHCHelperDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "\n",
    "eri_thc = kthc_eri_helper.get_eri(eri_kindices)\n",
    "# Can also do\n",
    "# eri_exact = kthc_eri_helper.get_eri_exact(eri_kindices)\n",
    "\n",
    "assert np.allclose(eri_thc, eri_exact)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's check convergence of the integral error, and corresponding MP2 error, with the THC dimension or equivalently the THC rank parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from openfermion.resource_estimates.pbc.hamiltonian.cc_extensions import compute_emp2_approx, build_cc_inst\n",
    "cc_inst = build_cc_inst(mf)\n",
    "eri_exact = cc_inst.ao2mo()\n",
    "emp2_exact, _, _ = cc_inst.init_amps(eri_exact)\n",
    "emp2_exact += mf.e_tot\n",
    "delta_eri = []\n",
    "delta_mp2 = []\n",
    "thc_ranks = np.arange(2, 20, 2)\n",
    "for cthc in thc_ranks:\n",
    "    num_thc = cthc * num_spatial_orbs\n",
    "    kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
    "    kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "    eri_thc = kthc_eri_helper.get_eri(eri_kindices)\n",
    "    eri_exact = kthc_eri_helper.get_eri_exact(eri_kindices)\n",
    "    # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
    "    delta_eri.append(np.max(np.abs(eri_thc-eri_exact))/num_kpts)\n",
    "    emp2_approx = compute_emp2_approx(mf, kthc_eri_helper)\n",
    "    delta_mp2.append(emp2_exact - emp2_approx)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the convergence of the integral error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(thc_ranks, delta_eri, marker=\"o\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$c_{\\mathrm{THC}}$\")\n",
    "plt.ylabel(r\"max$|\\Delta(pq|rs)|$\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how this corresponds to the MP2 error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.cla()\n",
    "plt.plot(thc_ranks, np.abs(delta_mp2), marker=\"o\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$c_{\\mathrm{THC}}$\")\n",
    "plt.ylabel(r\"$|\\Delta E_{\\mathrm{MP2}}|$ (Ha)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that apart from some non-monotonic behaviour (which is expected due to the non-linear nature of the ISDF procedure), that a relatively large rank parameter is required to obtain say $< 0.1$ mHa error per cell. Note this could likely be reduced by carefully selecting the orbital sets we perform ISDF on as oo, ov, and vv blocks exhibit different low-rank behaviour, but for quantum algorithms this is not relevant. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: Effect of Mesh Density\n",
    "\n",
    "You might be worried that we're cheating by only including 11^3 grid points, and that we're not saving much with the ranks we're choosing.\n",
    "\n",
    "Let us first see the fraction of points these ranks correspond to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(thc_ranks, 100*np.array(thc_ranks)*num_spatial_orbs/(11**3), marker=\"o\")\n",
    "plt.xlabel(r\"$c_{\\mathrm{THC}}$\")\n",
    "plt.ylabel(\"Percentage of real space points selected\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us crank up the FFTDF accuracy and see if the results change significantly. This cell will take around 10 minutes to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "results = {11: [], 15: [], 19: [], 21: [], 28: []}\n",
    "#results = {11: [], 15: []}\n",
    "for mesh in list(results.keys()):\n",
    "    ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
    "    cell = gto.Cell()\n",
    "    cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
    "    cell.a = ase_atom.cell[:].copy()\n",
    "    cell.basis = \"gth-szv\"\n",
    "    cell.pseudo = \"gth-hf-rev\"\n",
    "    # Use a smaller value of Ng that would otherwise be suggested (~ 26^3) for\n",
    "    # expediency + to allow exact ISDF factorization for comparison.\n",
    "    cell.mesh = [mesh]*3\n",
    "    cell.verbose = 0\n",
    "    cell.build()\n",
    "\n",
    "    kmesh = [1, 1, 3]\n",
    "    kpts = cell.make_kpts(kmesh)\n",
    "    num_kpts = len(kpts)\n",
    "    mf = scf.KRHF(cell, kpts)\n",
    "    mf.kernel()\n",
    "    print(\"SCF energy: \", mf.e_tot)\n",
    "\n",
    "    from pyscf.pbc.mp import KMP2 \n",
    "    emp2_exact = KMP2(mf).kernel()[0] + mf.e_tot\n",
    "    print(\"Ng = {}^3, MP2 Correlation energy: {}\".format(mesh, emp2_exact-mf.e_tot))\n",
    "    thc_ranks = np.arange(2, 20, 2)\n",
    "    for cthc in thc_ranks:\n",
    "        print(f\"Running mesh = {mesh}, cthc = {cthc}\")\n",
    "        num_thc = cthc * num_spatial_orbs\n",
    "        kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
    "        kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "        emp2_approx = compute_emp2_approx(mf, kthc_eri_helper)\n",
    "        # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
    "        results[mesh].append(emp2_exact - emp2_approx)\n",
    "\n",
    "    plt.plot(thc_ranks, np.abs(results[mesh]), marker=\"o\", label=f\"$N_g = {mesh}^3$\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$c_{\\mathrm{THC}}$\")\n",
    "plt.legend()\n",
    "plt.ylabel(r\"$|\\Delta E_{\\mathrm{MP2}}|$ (Ha)\")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Optional: Effect of Basis Set \n",
    "\n",
    "Another concern is that the basis set size is tiny so maybe things get worse as the basis set increases. Let's look into it by increasing the basis set size but still use a fairly coarse FFT grid, as we've seen it's not super important. This will also take several minutes to run.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "basis_results = {\"gth-szv\": [], \"gth-dzvp\": [], \"gth-tzvp\": []}\n",
    "for basis in list(basis_results.keys()):\n",
    "    ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
    "    cell = gto.Cell()\n",
    "    cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
    "    cell.a = ase_atom.cell[:].copy()\n",
    "    cell.basis = basis \n",
    "    cell.exp_to_discard = 0.1\n",
    "    cell.pseudo = \"gth-hf-rev\"\n",
    "    # Use a smaller value of Ng that would otherwise be suggested (~ 26^3) for\n",
    "    # expediency + to allow exact ISDF factorization for comparison.\n",
    "    cell.mesh = [15]*3\n",
    "    cell.verbose = 0\n",
    "    cell.build()\n",
    "\n",
    "    kmesh = [1, 1, 3]\n",
    "    kpts = cell.make_kpts(kmesh)\n",
    "    num_kpts = len(kpts)\n",
    "    mf = scf.KRHF(cell, kpts)\n",
    "    mf.kernel()\n",
    "    print(\"SCF energy: \", mf.e_tot)\n",
    "\n",
    "    num_spatial_orbs = mf.mo_coeff[0].shape[-1]\n",
    "    from pyscf.pbc.mp import KMP2 \n",
    "    emp2_exact = KMP2(mf).kernel()[0] + mf.e_tot\n",
    "    print(\"basis = {}, MP2 Correlation energy: {}\".format(basis, emp2_exact-mf.e_tot))\n",
    "    thc_ranks = np.arange(2, 20, 2)\n",
    "    for cthc in thc_ranks:\n",
    "        print(f\"Running basis = {basis}, cthc = {cthc}\")\n",
    "        num_thc = cthc * num_spatial_orbs\n",
    "        kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
    "        kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "        emp2_approx = compute_emp2_approx(mf, kthc_eri_helper)\n",
    "        # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
    "        basis_results[basis].append(emp2_exact - emp2_approx)\n",
    "\n",
    "    plt.plot(thc_ranks*num_spatial_orbs, np.abs(basis_results[basis]), marker=\"o\", label=f\"basis = {basis}\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$M$\")\n",
    "plt.legend()\n",
    "plt.ylabel(r\"$|\\Delta E_{\\mathrm{MP2}}|$ (Ha)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "for basis in list(basis_results.keys()):\n",
    "    plt.plot(thc_ranks, np.abs(basis_results[basis]), marker=\"o\", label=f\"basis = {basis}\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(r\"$c_{\\mathrm{THC}}$\")\n",
    "plt.legend()\n",
    "plt.ylabel(r\"$|\\Delta E_{\\mathrm{MP2}}|$ (Ha)\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effect on $\\lambda$ \n",
    "\n",
    "Now let us investigate the $\\lambda$ dependence of our ISDF-THC factorization. We will revert back to the minimal example from earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "ase_atom = bulk(\"C\", \"diamond\", a=3.5)\n",
    "cell = gto.Cell()\n",
    "cell.atom = pyscf_ase.ase_atoms_to_pyscf(ase_atom)\n",
    "cell.a = ase_atom.cell[:].copy()\n",
    "cell.basis = \"gth-dzv\"\n",
    "cell.pseudo = \"gth-hf-rev\"\n",
    "# Use a smaller value of Ng that would otherwise be suggested (~ 26^3) for\n",
    "# expediency + to allow exact ISDF factorization for comparison.\n",
    "cell.mesh = [15]*3\n",
    "cell.verbose = 0\n",
    "cell.build()\n",
    "\n",
    "kmesh = [1, 1, 3]\n",
    "kpts = cell.make_kpts(kmesh)\n",
    "num_kpts = len(kpts)\n",
    "mf = scf.KRHF(cell, kpts)\n",
    "mf.kernel()\n",
    "print(\"SCF energy: \", mf.e_tot)\n",
    "\n",
    "# converged SCF energy with appropriate Ng = -10.388904514046914, mesh = 28^3\n",
    "\n",
    "hcore = np.asarray([C.conj().T @ hc @ C for C, hc in zip(mf.mo_coeff, mf.get_hcore())])\n",
    "num_spatial_orbs = hcore.shape[-1]\n",
    "thc_ranks = np.arange(2, 20, 2)\n",
    "cc_inst = build_cc_inst(mf)\n",
    "eri_exact = cc_inst.ao2mo()\n",
    "emp2_exact, _, _ = cc_inst.init_amps(eri_exact)\n",
    "np.random.seed(7)\n",
    "print(f\"FFTDF MP2 Correlation energy: {emp2_exact}\")\n",
    "for cthc in thc_ranks:\n",
    "    num_thc = cthc * num_spatial_orbs\n",
    "    kpt_thc = kthc.solve_kmeans_kpisdf(mf, num_thc, verbose=False)\n",
    "    kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "    # Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
    "    thc_lambda = kthc.compute_lambda(hcore, kthc_eri_helper)\n",
    "    emp2_approx = compute_emp2_approx(mf, kthc_eri_helper) - mf.e_tot\n",
    "    emp2_error = abs(emp2_approx - emp2_exact)\n",
    "    print(f\"cthc = {cthc}, MP2 error: {emp2_error:4.3e}, lambda = {thc_lambda.lambda_total:4.3e}\")\n",
    "    if emp2_error < 1e-4:\n",
    "        print(f\"--> MP2 error < 0.1: {emp2_error:4.3e}, lambda = {thc_lambda.lambda_total:4.3e}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can try to improve $\\lambda$ by reoptimizing the THC factors using our ISDF factors as an initial guess. Practically this should mean we can use a smaller THC rank parameter for comparable MP2 accuracy. Note reoptimizing the THC factors is quite expensive. This cell may take 20 minutes to run. You should see that for a $c_\\mathrm{THC}=6$ the MP2 error is reduced by an order of magnitude. The optimization can be sped up by running on a GPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from openfermion.resource_estimates.pbc import utils\n",
    "# Recall we need RSGDF integrals to fit to.\n",
    "rsmf = scf.KRHF(mf.cell, mf.kpts).rs_density_fit()\n",
    "rsmf.kernel()\n",
    "mymp = KMP2(rsmf)\n",
    "emp2_rsgdf = mymp.kernel()[0]\n",
    "hcore, Luv = utils.build_hamiltonian(rsmf)\n",
    "np.random.seed(7)\n",
    "print(f\"RSGDF emp2: {emp2_rsgdf}\")\n",
    "cthc = 6\n",
    "num_thc = cthc * num_spatial_orbs\n",
    "# Here we use the helper function kpoint_thc_via_isdf which will first find\n",
    "# the ISDF factors and feed these into the BFGS and AdaGrad solvers.\n",
    "kpt_thc, _ = kthc.kpoint_thc_via_isdf(mf, Luv, num_thc, verbose=False)\n",
    "kthc_eri_helper = kthc.KPTHCDoubleTranslation(chi=kpt_thc.chi, zeta=kpt_thc.zeta, kmf=mf)\n",
    "# Note pyscf omits a normalization factor of 1/Nk in their definition of ERIs\n",
    "thc_lambda = kthc.compute_lambda(hcore, kthc_eri_helper)\n",
    "emp2_approx = utils.compute_emp2_approx(mf, kthc_eri_helper) - mf.e_tot # only compare correlation energy.\n",
    "emp2_error = abs(emp2_approx - emp2_rsgdf)\n",
    "print(f\"cthc = {cthc}, MP2 error: {emp2_error:4.3e}, lambda = {thc_lambda.lambda_total:4.3e}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyscf_pip",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "61e04ed291edad8fb55208ca4954976506ff083febf358c737688b49027371c1"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
